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(U)^This  report  is  concerned  with  the  application 
of  general  purpose  finite  element  structural  analysis 
programs  to  warhead  design  work.  The  capabilities  of 
a number  of  linear  and  nonlinear  programs  including 
NAS TRAN,  SAP  IV,  ADINA,  and  HONDO  are  described.  These 
programs  are  used  to  solve  a number  of  impact-type 
problems  ranging  in  complexity  from  a simple  one- 
dimensional rod  to  a conical  steel  warhead  filled  with 
explosive.  The  results  are  compared  on  the  basis  of 
accuracy,  time  and  cost,  ease  of  learning,  and  ease 
of  use. 
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INTRODUCTION 


Warheads  intended  to  penetrate  a target  and  detonate  within  it  must 
be  able  to  withstand  the  shock  loading,  resulting  from  impact.  In  partic- 
ular, the  warhead  case,  the  explosive,  and  the  fuze  must  each  remain  in- 
tact and  functional.  The  rational  design  of  such  warheads  requires  a 
knowledge  of  the  transient  stresses,  displacements,  and/or  accelerations 
produced.  These  can  then  be  compared  with  appropriate  damage  criteria 
to  determine  the  survivability  of  a particular  configuration.  Experi- 
mental data  on  these  variables  is  difficult  to  obtain.  Structural  anal- 
ysis programs  employing  the  finite  element  method  provide  the  designer 
with  an  attractive  alternative  method  of  obtaining  this  information. 

Although  the  finite  element  method  has  been  in  use  for  a number  of 
years,  it  is  only  comparatively  recently  that  codes  have  been  developed 
that  bring  the  solution  of  large-scale  transient  problems,  and  particu- 
larly nonlinear  problems,  to  the  threshold  of  economic  feasib ' 1 ity. 1 In 
the  finite  element  method,  the  body  to  be  analyzed,  which  may  possess  an 
infinite  number  of  degrees  of  freedom,  is  represented  by  an  approximate 
model  consisting  of  a number  of  idealized  structural  elements  connected 
at  specified  nodal  points,  and  possessing  only  a finite  and  relatively 
small  number  of  degrees  of  freedom.  The  response  of  this  model  to  the 
applied  loads  or  displacements  is  then  determined.  The  accuracy  of  the 
solution  obtained  depends  primarily  on  the  number  of  degrees  of  freedom 
in  the  model  and,  in  the  case  of  dynamic  problems,  on  the  time  step  used 
in  the  numerical  integration. 

The  limitations  of  structural  analysis  programs  for  design  work  are 
primarily  those  of  cost  and  time.  Impact  problems,  because  of  the  high 
stress  gradients  and  associated  steep  wave  fronts  involved,  require 
models  with  many  degrees  of  freedom  and  small  time  steps  to  obtain  accu- 
rate results.  The  cost  of  a finite  element  solution  can  be  expected  to 
increase  approximately  as  the  square  of  the  number  of  degrees  of  freedom 
and  directly  as  the  number  of  time  steps  used.  Not  only  does  the  com- 
puter cost  become  large  for  impact  problems,  but  the  time  in  man-hours 
for  the  construction  of  models  and  the  input  of  data  also  becomes  large. 


1 Finite  Element  Analysis  of  Transient  Nonlinear  Structural  Behavior , 
T.  Belytschko,  et  al,  eds.  Am er.  Soc.  Mech.  Eng.,  AMD-Vol.  14,  1975. 
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In  addition,  the  amount  of  output  data  generated  and  the  time  required 
for  its  interpretation  increase  significantly. 

% 

In  this  report  a number  of  linear  and  nonlinear  structural  analysis 
programs  are  described  and  then  applied  to  the  solution  of  impact-type 
problems  ranging  in  complexity  from  a simple  one-dlnenslonal  rod  to  a 
conical  warhead  filled  with  explosive.  The  results  provide  information 
as  to  the  suitability  of  finite  element  programs  for  warhead  design  work. 

The  various  programs  are  compared  on  the  basis  of  accuracy,  time  and 
cost,  ease  of  learning,  and  ease  of  use.  Recommendations  as  to  which 
programs  appear  most  useful  for  different  problems  are  offered. 


DESCRIPTION  OF  PROGRAMS 


A wide  variety  of  linear  and  nonlinear  structural  analysis  programs 
have  been  developed  and  are  being  used  in  government  and  Industry  at  the 
present  time.  Some  of  these  are  proprietary,  others  are  available  on  a 
lease  basis,  and  stiil  others  can  be  obtained  free  or  at  nominal  cost. 

The  particular  programs  looked  at  in  this  study,  most  of  which  fall  into 
the  last  category,  are  described  in  the  following  paragraphs.  In  addi- 
tion, some  features  of  the  various  programs  are  compared  in  Table  1. 

NASTRAN  is  probably  the  most  widely  used  and  continuously  developed 
program  available  for  linear  analysis.  Impact  problems  can  be  analyzed 
by  either  the  transient  or  modal  techniques.  A full  range  of  element 
types  is  available.  NASTRAN  has  excellent  routines  for  generating  struc- 
ture plots  and  plots  of  stresses  and  other  output  quantities.  However, 
a lack  of  grid  and  element  data  generation  routines  makes  problem  input 
time  consuming.  Initial  training  in  the  use  of  this  program  is  required 
because  of  the  large  number  of  options  available  and  a poor  method  of 
indexing  in  the  user's  manual. 

SAP  IV  is  another  widely  used  linear  analysis  program.  Its  range  of 
features  and  capabilities  is  not  quite  as  broad  as  in  the  case  of  NASTRAN. 
On  the  other  hand,  it  is  simpler  to  learn  and  easier  to  use,  especially 
for  the  type  of  problems  considered  in  this  study.  SAP  IV  has  convenient 
grid  and  element  data  generation  routines  that  reduce  the  amount  of  work 
involved  in  problem  input. 

SHELL  SHOCK  is  yet  another  linear  program  capable  of  handling  tran- 
sient problems.  It  offers  a choice  between  implicit  and  explicit  time 
integration  techniques.  Three  types  of  structural  elements  are  gener- 
ated internally:  beams,  shells,  and  solids  of  revolution.  Although  re- 
stricted to  axisymmetric  structures,  the  program  can  treat  nonaxisym- 
metric  loadings  (as  might  occur  during  the  oblique  impact  of  warheads). 
The  loading  is  broken  down  into  a series  of  Fourier  harmonics,  and  the 
response  to  each  harmonic  is  determined  separately. 
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NONSAP  is  a nonlinear  program  capable  of  solving  impact  problems  by 
either  modal  superposition  or  direct  Integration.  Both  geometric  and 
material  nonlinearities  (i.e.,  large  deformations  and  plasticity)  can 
be  handled.  However,  the  nonlinear  capabilities  are  limited  to  axisym- 
metric  configurations  subjected  to  axisymmerric  loads.  Like  SAP  IV,  to 
which  it  is  similar  in  format,  NONSAP  is  easy  to  use.  Currently,  the 
dynamic  capabilities  of  the  program  are  not  functional  on  the  NWC  UNIVAC 
1110,  so  only  static  analyses  can  be  run.  This  program  is  not  supported 
by  the  developers. 

ADINA  is  a follow-on  program  to  NONSAP  that  has  been  fully  imple- 
mented on  the  UNIVAC  machine.  ADINA  offers  additional  material  models, 
and  its  nonlinear  capabilities  extend  to  general  three-dimensional  prob- 
lems. 

HONDO  can  be  used  to  calculate  the  large  deformation  elastic,  or  in- 
elastic transient  response  of  two-dimensional  solids.  Time  integrations 
are  performed  using  an  explicit  central  difference  method.  Solutions 
are  obtained  entirely  in  core.  This  makes  the  program  very  fast  but  re- 
stricts the  model  size.  Eight  material  subroutines  are  available  includ- 
ing elastic,  elastic-plastic,  and  viscoelastic  behavior.  Though  somewhat 
limited  as  to  the  class  of  problems  it  can  solve,  HONDO  is  by  far  the 
simplest  of  the  programs  to  use. 

MARC  has  a reputation  for  the  broadest  capability  of  nonlinear  anal- 
ysis and  also  the  slowest  execution  time.  It  has  a full  complement  of 
modern  elements  and  material  representations.  It  can  treat  nonaxisym- 
metric,  nonlinear  transient  deformations.  There  are  excellent  pre  and 
post  processors  including  plotting  that  make  the  program  handy  to  use. 

It  seems  to  be  easy  to  learn.  MARC  is  a proprietary  program  obtained 
on  a lease  basis. 

TOODY-IIA  is  a two-dimensional  Lagrangian  hydrocode  suitable  for 
computing  wave  propagation  in  two  dimensions  in  rectangular  or  cylindri- 
cal coordinates.  The  code  is  based  on  the  conventional  f inite  difference 
analogs  to  the  nonlinear  Lagrangian  equations  of  motion.  TOODY-IIA  is 
applicable  over  a wide  range  of  conditions  including  such  material  pro- 
perties as  elastic-plastic,  crushable  foams,  deformable  solids,  explo- 
sives, and  fluids.  The  geometry  capability  is  very  general  and  col- 
liding objects  can  be  modeled.  This  code  provides  a sliding  interface 
option. 


TRANSIENT  LOADING  PROBLEMS 


ONE-DIMENSIONAL  ROD 

The  first  problem  considered  is  that  of  a one-dimensional,  elastic 
rod  built  in  at  one  end  and  loaded  at  the  other  end  with  a half  sine 


4 


NWC  TP  5981 


pressure  pulse.  The  rod  has  unit  properties,  while  the  pressure  pulse 
is  of  unit  maximum  intensity  and  is  applied  over  a period  of  time  such 
that  the  length  of  the  pulse  is  one-tenth  the  length  of  the  rod.  The 
problem  geometry  is  shown  in  Figure  1.  This  problem  is  included  because 
it  is  relatively  simple  and  has  a known  analytic  solution  and  not  be- 
cause it  is  a reasonable  representation  of  an  impacting  warhead. 

The  problem  was  solved  using  N'ASTRAN  and  SAr  IV.  For  modeling  pur- 
poses the  rod  was  divided  into  20  truss-type  elements  and  an  integration 
time  step  chosen  such  that  100  steps  were  required  for  the  pressure  pulse 
to  travel  the  length  of  the  rod  and  return  to  the  free  end.  Centerline 
displacements  along  the  length  of  the  rod  are  plotted  for  five  success- 
ive times  in  Figure  2.  NASTRAN  and  SAP  IV  results  are  shown  along  with 
the  exact  solution.  A zone  of  constant  displacement  is  created  behind 
the  pressure  pulse  as  it  travels  down  the  rod.  After  the  pulse  is  re- 
flected from  the  built-in  end,  this  displacement  is  removed  as  the  pulse 
travels  back  to  the  free  end.  Under  the  pressure  pulse  itself,  the  dis- 
placement varies  sinusoidally  from  zero  to  its  constant  value. 

Except  for  the  presence  of  oscillations  about  the  true  solution, 
both  programs  provide  a reasonably  good  gross  representation  of  the 
zones  of  zero  and  constant  displacement  in  the  rod.  On  the  other  hand, 
the  finite  element  solutions  do  not  accurately  describe  the  displacement 
under  the  pulse  itself.  Oscillations,  such  as  observed  here,  are  typi- 
cal of  finite  element  solutions  to  transient  problems,  and  in  more  com- 
plex situations  they  sometimes  render  the  interpretation  of  results 
difficult.  It  would  appear  that  the  amplitude  of  the  oscillations  is 
slightly  less  in  the  SAP  IV  solution  (20  elements),  possibly  due  to  a 
larger  artificial  viscosity  effect  from  the  Wilson  0 method  of  integra- 
tion used  in  that"  program. 


20  ELEMENTS 
0.02  TIME  STEP 


FIGURE  1.  One-Dimensional  Rod  Geometry. 
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The  accuracy  of  a finite  element  solution  can  be  improved  through 
the  use  of  more  elements  and  smaller  integration  time  steps.  Also  shown 
in  Figure  2 are  SAP  IV  results  obtained  using  a model  consisting  of  100 
elements  integrated  over  200  time  steps.  It  is  obvious  that  the  oscil- 
lations about  the  true  solution  in  the  zones  of  zero  and  constant  dis- 
placement have  been  substantially  reduced.  However,  the  accuracy  under 
the  pulse  itself  is  still  not  very  good,  especially  at  later  times. 

In  actual  design  work  the  shape  of  the  applied  pressure  pulse  would 
not  be  known  with  any  degree  of  precision,  and  the  exact  displacement 
under  the  pulse  would  not  be  of  interest.  Hence,  the  solution  given 
here  would  be  satisfactory.  If  greater  accuracy  were  desired,  however, 
a finer  model  consisting  of  about  1,000  elements  could  be  used.  For 
this  simple  one-dimensional  problem  the  additional  run  cost  would  be  in- 
significant. In  the  case  of  a two-  or  three-dimensional  problem,  how- 
ever, a model  with  the  equivalent  of  1,000  elements  in  two  or  three 
directions  would  be  extremely  expensive  in  terms  of  computer  time. 

A similar  problem  was  solved  for  a material  with  nonlinear  material 
properties.  In  particular,  bilinear  behavior  was  assumed  corresponding 
to  an  elastic-plastic  material  with  linear  strain  hardening.  The  moduli 
in  the  two  linear  regions  were  chosen  so  as  to  make  the  plastic  wave 
speed  half  the  elastic  wave  speed.  The  rod,  again,  had  unit  properties 
and  was  built-in  at  one  end.  In  this  case,  however,  the  other  end  was 
loaded  with  a step  pressure  equal  to  twice  the  yield  stress. 

This  problem  was  solved  using  ADINA  and  HONDO.  For  the  ADINA  solu- 
tion, the  rod  was  modeled  with  100  beam  elements.  A time  step  of  .01 
was  used.  Since  HONDO  does  not  have  one-dimensional  elements  in  its 
library,  it  was  necessary  to  use  a two-dimensional  representation  con- 
sisting of  a single  row  of  100  axisymmetric  elements  down  the  length  of 
the  rod.  To  minimize  transverse  inertia  effects,  the  diameter  of  the 
rod  was  made  small  compared  to  its  length  (D/L  = .01).  It  was  not  nec- 
essary in  this  case  to  select  an  integration  time  step  since  HONDO  auto- 
matically computes  an  appropriate  step  such  as  to  ensure  stability  of 
the  numerical  integration. 

The  results  are  shown  in  Figure  3,  which  is  a series  of  plots  of 
axial  stress  versus  axial  distance  from  the  free  end  at  selected  times. 
The  exact  solution,  which  is  characterized  by  the  presence  of  elastic 
and/or  plastic  wave  fronts  traveling  along  the  rod,  is  also  shown. 
Agreement  between  HONDO  and  the  exact  solution  is  quite  good  except  in 
the  immediate  vicinity  of  the  discontinuities  at  the  wave  fronts,  which 
any  finite  element  program  would  have  difficulty  representing.  The 
ADINA  solution,  however,  oscillates  wildly  around  the  exact  solution. 
These  oscillations  only  occur  in  the  region  of  plastic  behavior.  ADINA 
results  are  plotted  for  the  first  time  step  only  because  they  get  worse 
at  later  times.  It  is  not  clear  whether  these  oscillations  are  due  to 
an  error  in  the  version  of  the  program  implemented  at  NWC,  or  whether 
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they  are  inherent  in  the  program  itself,  perhaps  resulting  from  the  im- 
plicit integration  technique  used.  A new  version  of  ADINA  is  scheduled 
for  release  soon.  This  new  version  will  have,  among  other  features,  an 
explicit  time  integration  option.  It  would  be  of  interest  to  rerun  this 
problem  using  the  new  version  to  see  if  improved  results  are  obtained. 


ROD  WITH  TRANSVERSE  INERTIA 

In  the  one-dimensional  rod  idealization,  the  effects  of  transverse 
inertia  of  the  material  are  neglected!  Inclusion  of  these  effects  by 
means  of  the  Pochhammer-Chree  equations  greatly  complicates  the  problem. 
The  shape  of  a wave  front  propagating  down  the  rod  no  longer  remains 
constant,  but  becomes  degraded  and  dispersed  due  to  the  transfer  of  part 
of  its  energy  into  transverse  vibration.  Even  for  such  a simple  appear- 
ing problem  as  the  application  of  a step  pressure  to  the  end  of  a semi- 
infinite circular  cylinder,  no  satisfactory  analytic  solution  exists. 

A formal  solution  to  this  problem  as  a sum  of  Fourier  integrals  has  been 
obtained  by  Folk,  et  al.2  These  integrals  cannot  be  evaluated  by  simple 
means,  and  Folk,  et  al,  were  forced  to  obtain  asymptotic  solutions  which 
are  approximately  valid  only  at  large  distances  from  the  loaded  end. 

Noting  the  lack  of  analytical  results  for  cylinders  near  the  point 
of  application  of  a load.  Bertholf3  employed  a finite  difference  tech- 
nique to  solve  the  problem  of  a one-diameter-long  cylinder  subjected  to 
a step  pressure  at  one  end  and  rigid-lubricated  at  the  other.  Because 
his  solution  appears  to  be  quite  accurate,  the  same  problem  was  run  in 
this  study  using  NASTRAN  and  SAP  IV  and  the  results  compared.  A 300- 
element  model  of  the  cylinder  was  used.  The  cylinder  geometry  and  pro- 
perties, as  well  as  the  finite  element  representation,  are  shown  in 
Figure  4. 

Plots  of  axial  strain  at  the  free  surface  of  the  cylinder  versus 
axial  distance  from  the  loaded  end  for  lour  successive  nondimensional 
times  are  given  in  Figure  5.  The  Bertholf  solution,  the  two  finite  ele- 
ment solutions,  and  also  the  one-dimensional  rod  solution  are  shown.  It 
can  be  seen  that  the  two  finite  element  solutions  agree  quite  closely 
with  one  another  and  are  in  reasonably  good  overall  agreement  with  the 
Bertholf  solution.  However,  the  fine  points  of  the  Bertholf  solution 
(in  particular,  the  undulations  behind  the  front  of  the  wave)  are  not 


2 R.  Folk,  et  al.  "Elastic  Strain  Produced  by  Sudden  Application 
of  Pressure  to  One  End  of  a Cylindrical  Bar,"  Acoust.  Soa.  Amer. , J, 

Vol . 30  (1958),  p.  552. 

3 L.  D.  Bertholf.  "Numerical  Solution  for  Two-Dimensional  Elastic 
Wave  Propagation  in  Finite  Bars,"  J.  Appl.  Meoh. , Vol.  34,  Trans.,  ASME 
(Vol.  39),  Series  E (1967),  p.  725. 
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picked  up  by  either  finite  element  solution.  The  NASTRAN  solution  is 
in  slightly  better  agreement  with  the  Bertholf  solution,  especially  in 
its  representation  of  the  steepness  of  the  wave  front.  The  SAP  IV  solu- 
tion is  more  smoothed  out  and  less  detailed,  again  perhaps  due  to  a lar- 
ger damping  effect  from  the  Wilson  0 method  of  integration.  It  is  pos- 
sible that  the  Bertholf  solution,  being  itself  an  approximate  solution, 
is  no  more  accurate  or  correct  than  the  others.  Assuming  that  it  is, 
however,  it  follows  that  a considerably  more  refined  finite  element  model 
is  required  to  approximate  the  Bertholf  results.  It  should  be  noted  that 
the  one-dimensional  rod  solution  becomes  an  increasingly  poor  descrip- 
tion of  the  actual  strain  behavior  as'  time  increases. 

Bertholf  and  Karnes4  give  a finite  difference  solution  to  a similar 
problem,  this  time  involving  a nonlinear  material.  Specifically,  the 
problem  is  that  of  the  longitudinal  impact  of  two  elastic-plastic  cylin- 
ders. The  cylinders  were  0.0635  meters  long,  0.0318  meters  in  diameter, 
and  were  made  of  6061-T6  aluminum.  The  elastic-plastic  material  was 
assumed  to  follow  von  Mises'  yield  criterion  and  have  isotropic  strain 
hardening.  An  approach  velocity  of  the  two  cylinders  of  100  m/s  was 
considered.  Bertholf  and  Karnes  also  obtained  experimental  results  with 
which  to  compare  their  theoretical  solution. 

This  problem  was  solved  in  the  present  study  using  HONDO.  By  sym- 
metry, only  one  of  the  cylinders  was  considered.  A model  consisting  of 
six  layers  of  axisymmetric  elements  with  25  elements  in  each  layer  was 
used.  The  cylinder  was  given  an  initial  longitudinal  velocity  of  50  m/s. 
Points  at  the  impact  end  were  constrained  from  movement  in  the  longitud- 
inal direction,  but  were  allowed  to  move  freely  in  the  radial  direction 
(rigid-lubricated).  The  geometry  and  properties  are  shown  in  Figure  6. 

Results  in  the  form  of  plots  of  axial  strain  at  the  surface  versus 
time  at  distances  z = R/12,  R and  2R  from  the  impact  end  are  given  in 
Figure  7.  The  numerical  and  experimental  results  of  Bertholf  and  Karnes 
are  also  shown.  Agreement  between  theory  and  experiment  appears  to  be 
quite  good  considering  the  experimental  difficulties  involved.  The 
Bertholf  and  Karnes  solution  is  in  better  agreement  with  experiment  than 
the  HONDO  solution,  especially  at  larger  times.  This  is  to  be  expected 
since  the  Bertholf  and  Karnes  solution  utilizes  a much  finer  mesh  (3,600 
grid  points  as  opposed  to  182  for  HONDO).  However,  it  should  be  noted 
that  the  Bertholf  and  Karnes  solution  required  about  three  hours  on  a 
CDC  3600  computer  compared  to  about  one  minute  on  the  UNIVAC  1110  for 
the  HONDO  solution.  The  greater  error  in  the  HONDO  solution,  therefore, 
must  be  weighed  against  the  fact  that  it  cost  virtually  nothing  to  ob- 
tain. 


4 L.  D.  Bertholf  and  C.  H.  Karnes.  "Axisymmetric  Elastic-Plastic 
Wave  Propagation  in  6061-T6  Aluminum  Bars  of  Finite  Length,"  J.  Appl . 
Meah. , Vol.  36,  Trans.,  ASME,  Vol.  91,  Series  E (1969),  p.  533. 
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FIGURE  5.  (Contd 
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FIGURE  6.  Geometry  and  Properties  of  Bertholf-Karnes 
Elastic,  Plastic  Cylinder. 


EXPLOSIVE-FILLED  STEEL  CYLINDER 


Comparison  of  Codes.  Warhead  impact  with  a hard  target  at  normal 
obliquity  has  been  investigated  with  the  model  shown  in  Figure  8.  A 0.1- 
meter-long  cylinder,  0.1  meter  in  diameter  is  composed  of  a thin  steel 
case  filled  with  explosive.  The  case  is  assumed  to  be  elastic  for  stress 
states  within  the  von  Mises  yield  condition  and  have  linear  kinematic 
strain  hardening  for  stress  states  satisfying  the  yield  criterion.  The 
explosive  is  assumed  to  be  a linear  elastic  material.  Material  properties 
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FIGURE  7.  Elastic-Plastic  Rod  With 
Transverse  Inertia,  Axial  Strain  at 
Free  Surface  Versus  Axial  Distance. 
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used  in  this  analysis  are  shown  in  Table  2.  The  failure  stress  for  steel 
is  the  one-dimensional  tensile  yield  stress  while  the  explosive  failure 
stress  is  the  lowest  pressure  at  which  deflagration  occurs.  A coarse 
finite  element  model  of  this  warhead  configuration  has  been  analyzed. 

The  model  uses  three  layers  of  axisymraetric  elements  for  the  explosive 
core  and  one  layer  for  the  case  (see  Figure  8) . 


TABLE  2.  Cylinder  Material  Properties. 


Material 

Stress-strain 

constants 

Density 
p,  Kg/m3 

Wave  speed 
m/sec 

Failure- 

E,  Pa 

U[ | 

stress , 
Pa 

Steel 

(elastic) 

2.0  x 10n 

0.3 

7800 

Cr  = 5060 

6 . 9 x 108 

Steel 

(plastic) 

2.0  x 109 

0.3 

7800 

Cr  = 506 

— 

PBX 

(expl) 5 

6.9  x 109 

0.34 

1800 

Ci  = 2420 

5. Ox  108 

A half  sine  pressure  pulse  of  5 usee  duration  is  applied  to  one  end 
of  the  cylinder.  Two  loading  cases  wii_h  different  pressure  amplitudes 
have  been  defined  in  Table  3.  Load  1 should  impart  a uniform  1 m/sec 
axial  velocity  to  the  end  of  the  cylinder.  The  resultant  stresses  at  the 
impacting  end  will  be  in  the  plastic  range.  Only  those  codes  with  the 
capability  to  analyze  nonlinear  material  response  (ADINA,  HONDO  and  MARC) 
have  been  used  for-  case  2. 


TABLE  3.  Pressure  Amplitude. 


Case 

P]^  (explosive) 

P2  (case) 

1 

4.32  x 106 *  Pa 

10.8  x 10s  Pa 

2 

39.46  x 106  Pa 

16.06  x 108  Pa 

Axial  displacements  on  the  centerline  and  on  the  explosive-case  in- 
terface have  been  calculated  by  the  appropriate  codes.  These  displace- 
ments are  compared  at  10,  20,  and  30  psec  after  impact  in  Figures  9 and 
10  for  elastic  and  elastic-plastic  deformations,  respectively. 


5 Lawrence  Livermore  Laboratory.  Properties  of  Chemical  Explosives 

and  Explosive  Simulants , by  B.  M.  Dobratz.  University  of  California, 

July  31,  1974.  (UCRL  51319,  Rev.  1.) 
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FIGURE  9.  Comparison  of  Linear  Analyses  for  80-Element  Cylinder  Impact 
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Within  each  load  case  all  of  the  programs  predict  similar  axial  dis- 
placements on  the  case-explosive  interface.  Agreement  on  the  centerline 
displacements  is  not  as  good  and  the  differences  between  the  solutions 
increase  with  increasing  time. 

Cylinder  Impact.  Rather  than  the  half  sine  pressure  pulse  that  was 
used  to  compare  codes,  impact  loading  can  be  represented  as  an  initial 
velocity  condition.  Normal  end-on  impact  with  a rigid  wall  requires  a 
uniform  axial  initial  velocity  of  the  cylinder  together  with  zero  dis- 
placement at  the  impacting  end.  NASTRAN  was  used  to  obtain  the  early 
time  response  of  the  model  described  in  Figure  8.  Both  stresses  and 
displacements  were  output  for  times  up  to  80  psec  after  impact.  Typi- 
cal axial  and  circumferential  strain  histories  on  the  case-explosive 
interface,  and  on  the  centerline  are  shown  in  Figures  11  through  13, 
respectively. 

In  Figure  11  a compressive  wave  is  observed  to  be  traveling  along 
the  length  of  the  cylinder  at  approximately  5 mm/psec.  It  takes  about 
20  psec  for  the  wave  to  traverse  the  100-mm-long  cylinder.  A reflected 
wave  from  the  free  end  of  the  cylinder  does  not  completely  relieve  the 
strains  because  of  interactions  between  the  case  and  the  explosive  core. 
Circumferential  strains  in  the  case  (see  Figure  12)  are  caused  by  both 
the  Poisson  effect  and  interface  pressures.  In  this  relatively  thin 
case,  interface  pressures  predominate.  Consequently,  the  circumferen- 
tial strains  travel  at  a slower  speed,  corresponding  to  the  wave  propa- 
gation speed  in  the  less  rigid  explosive  material. 

The  strain  history  on  the  centerline  of  the  explosive  is  less  smooth 
than  that  on  the  case  (Figure  13).  The  explosive,  being  less  rigid  than 
the  case,  is  similar  to  a fluid.  In  addition  to  an  irrotational  wave, 
it  responds  to  the  wall  motion.  The  apparent  wave  speed  of  3 mm/psec  is 
somewhat  faster  than  the  irrotational  wave  speed  (2.4  mm/psec). 

Calculated  strains  can  be  compared  with  predictions  based  on  a one- 
dimensional wave  propagation  theory.  In  the  explosive,  plane  strain 
should  predominate.  The  corresponding  stress,  based  on  quantities  from 
Table  2 is 

ozz/V  = pCi  = 4.3  x 106  Pa  sec/m 

where  V is  the  impact  velocity.  For  the  steel  case,  a plane  stress 
solution  yields 

oz2/V  = pCr  = 39.4  x 106  Pa  sec/m 
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FIGURE  11.  Cylinder  Impact  With  Rigid  Wall  (At  * 0.5  x 10~6  sec,  r/R 


FIGURE  12.  Cylinder  Impact  With  Rigid  Wall  (At  = 0.5  x 10-6  sec,  r/R  = 1/10). 
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AXIAL  COORDINATE,  m 


FIGURE  13.  Impact  With  Rigid  Wall  at  1 m/sec 
(At  = 0.5  x 10-6  sec,  r/R  = 1/10). 
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In  thin  cylinders,  the  hoop  stress  caused  by  hydrostatic  pressure  from 
the  explosive  would  be 

oqq/V  = pCt  R/T  = 43  x 106  Pa  sec/m  for  R/T  = 10 

Where  R is  the  cylinder  diameter  and  T is  the  case  thickness.  Corres- 
ponding strains  are  shown  in  Figures  11-13.  Generally,  maximum  stresses 
are  close  to  the  plane  wave  values  and  occur  near  the  impacted  end. 

Failure  Analysis.  Normal  impact  of  an  explosive-filled  cylinder 
with  a rigid  wall  can  result  in  failuire  by  breakup  of  the  case,  by  ex- 
plosive deflagration,  or  by  detonation  of  the  explosive.  Case  breakup 
and  deflagration  are  stress-dependent  failure  modes  while  detonation  is 
dependent  on  stress  history  in  the  explosive.  An  analysis  of  stresses 
in  elastic  cylinders  has  been  conducted  to  determine  how  the  minimum  im- 
pact speed  for  failure  depends  on  the  wall  thickness.  Surprisingly,  the 
lowest  predicted  speed  at  which  yield  occurs  in  a steel  case  and  deflag- 
ration occurs  in  the  explosive  core  are  both  only  weakly  dependent  on 
case  thickness. 

The  finite  element  analysis  described  in  the  previous  section  pro- 
vides stresses  in  each  element  used  to  represent  the  body.  From  an 
average  taken  from  a cross  section  near  the  impacted  end,  impact  stresses 
are  obtained  for  comparison  with  the  yield  and  deflagration  failure  cri- 
terion. Detonation  of  the  explosive  is  not  considered  since  we  are 
looking  for  the  lowest  impact  speed  at  which  failure  occurs.*  Both 
plastic  yielding  of  the  steel  cylinder  and  deflagration  of  the  explosive 
are  considered.  A 5-Kbar  pressure  deflagration  condition  has  been  as- 
sumed. The  resulting  minimum  impact  velocities  for  failure  are  shown  in 
Figure  14.  These-  failure  speeds  are  surprisingly  insensitive  to  the 
wall  thickness.  There  is  also  relatively  little  effect  of  the  cylinder 
end  conditions  on  the  failure  speeds.  End  conditions  considered  were 
that  the  impact  end  of  the  cylinder  was  either  free  to  expand  radially 
or  was  fixed  radially.  These  cases  simulate  impact  of  open  and  closed 
ended  warheads,  respectively. 

CONICAL  WARHEAD 

One  advantage  of  the  finite  element  procedure  is  its  ability  to 
handle  non-simple  geometries.  To  illustrate  this  capability,  another 


An  energy  density  failure  criterion  is  appropriate  for  detonation. 
An  example  of  a problem  for  application  of  this  criterion  would  be  de- 
termining the  maximum  plate  thickness  that  could  be  penetrated  at  some 
impact  speed  without  warhead  detonation. 
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warhead  configuration  consisting  of  an  explosive-filled  conical  steel 
case  with  linearly  varying  thickness  was  considered.  Two  different 
material  characterizations  were  used.  In  the  first,  both  the  case  and 
explosive  were  assumed  to  behave  elastically.  In  the  second,  the  ex- 
plosive was  assumed  to  behave  elastically,  but  the  case  was  treated  as 
an  elastic,  perfectly  plastic  material.  As  in  the  concentric  cylinder 
problem,  the  loading  consisted  of  a half  sine  pressure  pulse  at  one  end. 
The  other  end  was  built  in. 

Several  different  model  representations  were  used  in  the  runs  made 
with  the  various  programs.  For  the  NASTRAN  and  SAP  IV  runs,  four  layers 
of  axisymmetric  elements  were  used;  two  to  represent  the  case,  and  two 
to  represent  the  explosive.  There  were  20  elements  down  the  length  of 
the  warhead  in  each  layer.  The  integration  step  size  was  selected  as 
one-twentieth  of  the  pulse  length.  The  geometry  and  properties  of  this 
configuration  are  shown  in  Figure  15.  In  the  ADINA  runs,  six  layers 
(three  for  the  case,  and  three  for  the  explosive)  were  used.  The  HONDO, 
SHELL  SHOCK  and  TOODY-IIA  runs  were  made  with  variable  sized  meshes  with 
many  small  elements  near  the  loaded  end  and  fewer  relatively  large  ele- 
ments at  the  other  end.* 

The  results  are  summarized  in  Figures  16-19,  which  are  plots,  re- 
spectively, of  centerline  displacement,  interface  axial  displacement, 
case  axial  stress  near  the  interface,  and  case  radial  stress  near  the 
interface  versus  axial  distance  from  the  loaded  end  at  10  psec.  Sepa- 
rate plots  for  the  linear  and  nonlinear  runs  are  given  in  each  figure. 

It  is  obvious  that  there  are  considerable  differences  in  the  results  for 
the  various  programs,  due  in  part  to  differences  in  the  programs  them- 
selves (different  elements,  integration  schemes,  etc.),  and  also  to 
differences  in  the  model  representations  used.  Nevertheless,  the  gross 
behavior  of  the  various  solutions  is  reasonably  similar.  This  is  espe- 
cially true  for  the  displacements.  There  is  more  variation  in  the 
stress  results.  As  would  be  expected,  the  nonlinear  plots  show  larger 
peak  displacements  and  smaller  maximum  stresses  than  the  linear  ones. 

Although  the  results  for  the  various  programs  at  10  psec  are  in 
reasonable  agreement  with  each  other  (and  probably  with  the  "true"  solu- 
tion, if  one  were  available),  agreement  at  later  times  is  not  as  good. 
Figure  20  is  a series  of  plots  of  centerline  displacement  versus  axial 
distance  from  the  loaded  end  at  10  psec  intervals  up  to  80  psec  for  the 
NASTRAN  and  SAP  IV  runs.  After  about  40  psec,  it  is  difficult  to  recog- 
nize the  two  solutions  as  applying  to  the  same  problem.  This  indicates 


* The  solutions  to  this  problem  using  SHELL  SHOCK,  HONDO,  and  TOODY- 
IIA  were  obtained  by  Kaman  Sciences  Corp.  SHELL  SHOCK  and  TOODY-IIA  are 
not  available  on  the  NWC  UNIVAC  1110  at  the  present  time. 
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FIGURE  15.  Geometry  and  Properties  for  Conical  Warhead  Problem. 
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FIGURE  16.  Conical  Warhead  Impact,  Interface  Axial  Displacement  Versus  Axial  Distance 
Po  = 3.212  x 109  Pa,  t = 10  psec. 


FIGURE  17.  Conical  Warhead  Impact,  Interface  Axial  Displacement  Versus  Axial 
Distance,  p„  = 3.212  x 10  9 Pa,  t = 10  ysec. 
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that  the  model  used  is  not  sufficiently  detailed  to  obtain  accurate 
long-term  solutions.  A finer  grid  and  small  time  steps  would  be  re- 
quired to  produce  such  solutions. 

Breakdown  of  the  conical  warhead  solution  is  more  rapid  than  solu- 
tion breakdown  in  the  problems  previously  treated.  This  is  probably 
due  to  the  fact  that  in  the  previous  problems  the  primary  effect  of  the 
pressure  load  is  to  produce  longitudinal  stress  waves  which  travel  down 
the  length  of  the  cylinder.  Transverse  dispersion  effects  are  of  only 
secondary  importance.  In  the  conical  warhead  problem,  however,  due  to 
the  greater  case  thickness  and  conical  geometry,  transverse  waves  pro- 
duced by  reflection  from  the  conical  surfaces  will  make  an  important 
contribution  right  from  the  start  and  can  be  expected  to  cause  differ- 
ences between  the  solutions  at  earlier  times. 


SPHERICAL  CAVITY 

There  are  surprisingly  few  dynamic  plasticity  problems  involving 
more  than  a one-dimensional  state  of  stress  for  which  analytic  solutions 
are  available.  Even  such  a simple-appearing  problem  as  that  of  an  in- 
finite elastic,  perfectly  plastic  medium  containing  a spherical  cavity 
subjected  to  a step  pressure  cannot  be  solved  completely  by  analytic 
means.  However,  an  analytic  solution  to  this  problem  can  be  obtained 
for  short  times,  and  a good  numerical  solution  can  be  generated  for  lar- 
ger times.  It  was  decided,  therefore,  to  use  this  problem  to  test  the 
relative  efficiency  and  accuracy  of  some  of  the  nonlinear  structural 
analysis  programs. 

The  spherical  cavity  problem  for  an  elastic,  perfectly  plastic  me- 
dium was  discussed  at  length  in  the  survey  article  by  Hopkins.6  He 
gives  a formal  solution  in  terms  of  a set  of  ordinary  differential  equa- 
tions involving  the  plastic  wave  functions,  H and  K,  the  elastic  wave 
function,  F,  and  the  speed  of  the  elastic-plastic  boundary,  6. 

Application  of  a step  pressure  at  the  cavity  radius  causes  an  elas- 
tic wave  to  propagate  outward  into  the  medium  with  a constant  speed,  ce. 
If  the  pressure  is  greater  than  that  required  to  initiate  yielding,  then 
the  elastic  wave  will  be  followed  by  a plastic  wave  traveling  at  a slower 
but  constant  speed,  Id  = Cp.  There  will  be  a discontinuity  in  stress 
across  the  boundary  between  the  elastic  and  plastic  regions.  As  time 
increases,  the  magnitude  of  this  discontinuity  will  decrease  until  at 


6 H.  G.  Hopkins.  Progress  in  Solid  Mechanics,  North-Holland  Pub- 
lishing Co.,  Amsterdam,  Vol.  1,  Chapter  III,  1960,  pp.  83-164. 
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FIGURE  20.  Conical  Warhead  Problem.  Plots  of  centerline  displacement 
versus  axial  distance  at  selected  times  for  NASTRAN  and  SAP  IV  runs. 
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some  time,  tm,  it  will  disappear  entirely.  Up  to  this  time  an  analytic 
solution  to  the  problem  can  be  obtained  by  integration  of  the  governing 
equations  given  by  Hopkins. 

For  times  greater  than  tm,  the  character  of  the  solution  changes 
somewhat.  The  speed  of  the  elastic-plastic  boundary  no  longer  remains 
constant,  but  begins  to  decrease.  Because  the  boundary  speed  is  no 
longer  a known  constant,  an  additional  unknown  is  introduced  into  the 
problem.  Moreover,  because  this  unknown  enters  into  the  governing  equa- 
tions in  a nonlinear  way,  an  analytic  solution  is  no  longer  possible. 
However,  the  governing  equations  can  still  be  integrated  numerically 
using  an  iterative  procedure  to  eliminate  the  stress  discontinuity  be- 
tween the  elastic  and  plastic  regions.  Though  only  approximate,  this 
iterative  solution  can  be  expected  to  be  considerably  more  accurate 
than  finite  element  solutions  because  it  takes  into  account  the  possible 
presence  of  discontinuities  at  the  elastic  wave  front  and  the  elastic- 
plastic  boundary. 

This  solution  was  used  as  a standard  of  comparison  for  finite  ele- 
ment results  obtained  using  ADINA  and  HONDO.  The  spherical  cavity 
problem  is  one-dimensional.  Unfortunately,  none  of  the  structural 
analysis  programs  considered  here  have  one-dimensional  spherical  elements 
in  their  libraries.  Consequently,  it  was  necessary  to  use  a two-dimen- 
sional finite  element  model.  A narrow  conical  sliver  of  material  with 
the  apex  of  the  cone  coinciding  with  the  center  of  the  cavity  was  selec- 
ted for  analysis.  This  conical  sliver  was  modeled  using  axisymraetric 
elements.  Because  of  symmetry,  points  on  the  side  of  the  cone  can  only 
move  radially  outward.  This  boundary  condition  cannot  be  imposed  dir- 
ectly using  program  ADINA.  It  was  necessary,  therefore,  in  the  ADINA 
runs  to  attach  very  long,  stiff  and  massless  truss  elements  perpendic- 
ular to  the  side  at  these  points  to  prevent  radial  motion.  The  model 
geometry  is  shown  in  Figure  21. 

A comparison  of  results  obtained  using  ADINA  and  HONDO  with  the 
analytic-numeric  solution  is  given  in  Figure  22,  which  is  a series  of 
plots  of  radial  stress  versus  radial  distance  at  selected  times.  Agree- 
ment between  the  HONDO  results  and  the  "exact"  solution  is  quite  good, 
even  in  the  vicinity  of  the  stress  discontinuities,  which  the  program 
does  its  best  to  represent.  The  ADINA  results  are  not  as  satisfactory. 
Not  only  are  there  greater  errors  near  the  discontinuities,  but  fairly 
substantial  oscillations  about  the  true  solution  are  present.  Hope- 
fully, the  new  version  of  ADINA  due  to  be  released  soon  will  provide 
better  results  for  this  problem. 
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FIGURE  21.  Geometry  of  Spherical  Cavity  Model. 


CONCLUSIONS 


Based  on  the  foregoing  results,  the  following  conclusions  can  be 
drawn  regarding  the  usefulness  of  finite  element  structural  analysis 
programs  for  warhead  design  work. 


ACCURACY 

For  simple  geometries  where  the  wave  propagation  is  essentially  lon- 
gitudinal (e.g.,  the  rod  problems  and  the  Bertholf  cylinder)  reasonably 
accurate  solutions  were  obtained  with  the  mesh  sizes  and  time  steps 
used.  The  results  remained  adequate  for  the  time  required  for  a wave 
to  travel  the  length  of  the  cylinder  and  return  to  the  loaded  end.  In 
the  case  of  the  more  geometrically  complicated  conical  warhead,  the  re- 
sults became  essentially  meaningless  before  the  wave  had  traveled  even 
the  length  of  the  warhead.  Successful  analysis  of  this  problem  for  a 
longer  period  of  time  would  require  a much  finer  model  and  smaller  time 
steps.  The  cost  of  such  an  analysis  might  become  extremely  high. 

It  is  difficult  to  assess  the  relative  accuracy  of  the  various  pro- 
grams in  quantitative  terms  because  of  a lack  of  exact  solutions  (except 
for  the  very  simple  problems)  with  which  to  compare  the  results.  Clearly, 
the  version  of  ADINA  implemented  at  NWC  was  not  as  accurate  as  HONDO  in 
the  problems  to  which  they  were  both  applied.  On  the  other  hand,  while 
there  were  differences  in  the  results  obtained  using  the  other  programs, 
no  one  program  appeared  to  be  significantly  more  or  less  accurate  than 
the  others. 
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FIGURE  22.  Spherical  Cavity  Problem.  Radial  stress  versus  radial  distance 
at  selected  times. 
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TIME  AND  COST 

The  CAU  times  and  costs  on  the  NWC  UNIVAC  1110  system  for  the  vari- 
ous problem  solutions  are  given  in  Table  4.  For  some  of  the  runs  not 
made  at  NWC,  times  were  converted  to  equivalent  UNIVAC  1110  times.  The 
table  shows  that,  by  and  large,  those  programs  with  more  features  and 
capabilities  take  longer  to  solve  problems.  This  can  be  seen  by  com- 
paring the  NASTRAN  and  SAP  IV  solution  times  for  those  problems  they 
were  both  used  to  solve.  The  Bertholf  cylinder  problem,  for  example, 

required  596  CAU  seconds  and  cost  $255  using  NASTRAN compared  to  only 

182  seconds  and  $120  using  SAP  IV.  Because  of  its  smaller  size  (number 
of  FORTRAN  statements) , SAP  IV  was  able  to  solve  some  problems  in-core 
without  the  use  of  external  storage.  It  is  likely  that  for  extremely 
large  problems,  which  neither  program  could  solve  in-core,  NASTRAN  would 
come  into  its  own,  so  to  speak,  and  the  cost  and  time  differences  be- 
tween the  two  would  largely  disappear.  Surprisingly,  MARC,  which  has  a 
reputation  for  being  slow,  actually  required  less  CAU  time  than  did 
NASTRAN  to  solve  the  concentric  cylinder  problem.  HONDO,  which  is  an 
in-core  solver,  was  by  far  the  fastest  of  the  nonlinear  programs.  For 
Lhe  spherical  cavity  problem,  HONDO  required  only  55  CAU  seconds  and 
$14 compared  to  557  seconds  and  $158  for  ADINA. 


EASE  OF  USE 

For  the  most  part,  those  programs  requiring  the  least  amount  of  in- 
put tend  to  be  the  easiest  to  use.  Less  input  means  less  time  spent  on 
problem  preparation  and  less  likelihood  of  mistakes.  From  this  stand- 
point, it  would  appear  that  SAP  IV  should  be  the  program  of  choice  for 
three-dimensional  linear  elastic  problems  and  that  HONDO  should  be  used 
for  linear  and  nonlinear  problems  involving  axisymmetric  solids.  Use 
of  other  programs  is  justified  only  when  their  special  capabilities 
(such  as  the  fluid  elements  available  in  NASTRAN)  permit  the  solution 
of  problems  not  tractable  with  the  other  programs. 


EASE  OF  LEARNING 

The  difficulty  involved  in  learning  to  use  a program  depends  in 
large  measure  on  the  complexity  of  the  program  itself  and  on  the  amount 
of  input  required.  Hence,  programs  that  are  easy  to  use  tend  to  be  easy 
to  learn  how  to  use.  HONDO,  NONSAP  and  ADINA  are  offshoots  of  SAP  IV. 

A designer  or  analyst  familiar  with  the  format  and  nomenclature  used  in 
SAP  IV  can  expect  to  learn  to  use  these  other  programs  relatively 
quickly.  NASTRAN  is  especially  difficult  to  learn  because  the  user's 
manual  is  cumbersome  and  hard  to  read.  However,  once  the  manual  has 
been  mastered,  use  of  the  program  is  not  difficult. 
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TABLE  4.  Cost  and  Time  Comparison  for  the  Various  Runs. 


Problem 

Degrees  of 
freedom 

Bandwidth 

Number  of 
time  steps 

Program 

CAU  time 
seconds 

Cost , 
dollars 

1-D  rod,  elastic,  half 

20 

’2 

100 

NASTRAN 

19 

sine  pressure  pulse 

SAP  IV 

2 

| | 

100 

2 

200 

SAP  IV 

7 

7 

1-D  rod,  elastic-plastic. 

400 

8 

236 

HONDO 

57 

14 

step  pressure 

Bertholf  cylinder,  elastic. 

640 

25 

80 

NASTRAN 

596 

step  pressure 

SAP  IV 

182 

Bertholf  cylinder,  elastic- 

364 

HONDO 

33 

12 

plastic,  impact 

Concentric  cylindrical  war- 

110 

14 

30 

NASTRAN 

60 

head,  elastic,  half  sine 

NONSAP  (In) 

3 

pressure  pulse 

210 

14 

30 

NASTRAN 

W 

SAP  IV 

mm 

NONSAP  (In) 

WtM 

546 

30 

30 

SAP  IV 

44 

Concentric  cylindrical  war- 

110 

14 

30 

NONSAP 

41 

head,  elastic-plastic. 

210 

14 

30 

NONSAP 

75 

half  sine  pressure  pulse 

MARC 

175 

ADINA 

150 

116 

HONDO 

8 

3 

Conical  warhead,  elastic 

176 

13 

80 

NASTRAN 

El 

mm 

half  sine  pressure  pulse 

SAP  IV 

El 

■a 

- 

510 

21 

160 

SAP  IV 

223 

sa 

Conical  warhead,  elastic. 

176 

13 

20 

NASTRAN 

95 

52 

1 m/s  impact 

Spherical  cavity,  elastic. 

302 

6 

88 

SAP  IV 

El 

step  pressure 

ADINA  (In) 

404 

8 

170 

HONDO  (In) 

30 

Spherical  cavity,  elastic- 

452 

6 

100 

ADINA 

557 

158 

plastic,  step  pressure 

604 

8 

150 

HONDO 

55 

14 

40 


NWC  TP  5981 


SUMMARY 


This  report  has  been  concerned  with  the  application  of  general  pur- 
pose finite  element  structural  analysis  programs  to  warhead  design  work. 
A number  of  linear  and  nonlinear  codes  were  described.  These  were 
applied  to  the  solution  of  impact-type  problems  of  varying  complexity. 
Based  on  the  results  obtained,  the  various  programs  were  compared  with 
regard  to  accuracy,  time  and  cost,  ease  of  learning,  and  ease  of  use. 
Considering  only  those  programs  currently  available  at  NWC,  it  would 
appear  that  SAP  IV  is  the  preferred  program  for  linear  elastic  problems, 
while  HONDO  is  best  suited  for  nonlinear  problems  involving  axial  sym- 
metry. Because  of  accuracy  problems,  the  current  version  of  ADINA 
should  be  used  with  caution. 
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